The present invention relates to a method and apparatus which uses model-based adaptive filtering techniques to estimate physiological parameters. More specifically, the invention employs Kalman filtering techniques in pulse oximetry to estimate the oxygen saturation of hemoglobin in arterial blood.
Pulse oximeters typically measure and display various blood flow characteristics including but not limited to the oxygen saturation of hemoglobin in arterial blood. Oximeters pass light through blood perfused tissue such as a finger or an ear, and photoelectrically sense the absorption of light in the tissue. The amount of light absorbed is then used to calculate the amount of the blood constituent (e.g., oxyhemoglobin) being measured.
The light passed through the tissue is selected to be of one or more wavelengths that are absorbed by the blood in an amount representative of the amount of the blood constituent present in the blood. The amount of light passed through the tissue varies in accordance with the changing amount of blood constituent in the tissue and the related light absorption.
When the measured blood parameter is the oxygen saturation of hemoglobin, a convenient starting point assumes a saturation calculation based on Lambert-Beer""s law. The following notation will be used herein:
I(xcex,t)=Io(xcex)exp(xe2x88x92(sxcex2o(xcex)+(1xe2x88x92s)xcex2r(xcex))1(t))xe2x80x83xe2x80x83(1)
where:
xcex=wavelength;
t=time;
I=intensity of light detected;
Io=intensity of light transmitted;
s=oxygen saturation;
xcex2o, xcex2r=empirically derived absorption coefficients; and
l(t)=a combination of concentration and path length from emitter to detector as a function of time.
The traditional approach measures light absorption at two wavelengths, e.g., red and infrared (IR), and then calculates saturation by solving for the xe2x80x9cratio of ratiosxe2x80x9d as follows.
1. First, the natural logarithm of (1) is taken (xe2x80x9clogxe2x80x9d will be used to represent the natural logarithm) for IR and Red
logI=logIoxe2x88x92(sxcex2o+(1xe2x88x92s)xcex2r)1xe2x80x83xe2x80x83(2)
2. (2) is then differentiated with respect to time                                                         ⅆ              log                        ⁢                          xe2x80x83                        ⁢            I                                ⅆ            t                          =                              -                          (                                                s                  ⁢                                      xe2x80x83                                    ⁢                                      β                    o                                                  +                                                      (                                          1                      -                      s                                        )                                    ⁢                                      β                    r                                                              )                                ⁢                                    ⅆ              l                                      ⅆ              t                                                          (        3        )            
3. Red (3) is divided by IR (3)                                                         ⅆ              log                        ⁢                          xe2x80x83                        ⁢                                          I                ⁡                                  (                                      λ                    R                                    )                                            /                              ⅆ                t                                                                        ⅆ              log                        ⁢                          xe2x80x83                        ⁢                                          I                ⁡                                  (                                      λ                    IR                                    )                                            /                              ⅆ                t                                                    =                                            s              ⁢                              xe2x80x83                            ⁢                                                β                  o                                ⁡                                  (                                      λ                    R                                    )                                                      +                                          (                                  1                  -                  s                                )                            ⁢                                                β                  r                                ⁡                                  (                                      λ                    R                                    )                                                                                        s              ⁢                              xe2x80x83                            ⁢                                                β                  o                                ⁡                                  (                                      λ                    IR                                    )                                                      +                                          (                                  1                  -                  s                                )                            ⁢                                                β                  r                                ⁡                                  (                                      λ                    IR                                    )                                                                                        (        4        )            
4. Solving for s   s  =                                                        ⅆ              log                        ⁢                          xe2x80x83                        ⁢                          I              ⁢                              (                                  λ                  IR                                )                                                          ⅆ            t                          ⁢                              β            r                    ⁢                      (                          λ              R                        )                              -                                                  ⅆ              log                        ⁢                          xe2x80x83                        ⁢                          I              ⁢                              (                                  λ                  R                                )                                                          ⅆ            t                          ⁢                              β            r                    ⁢                      (                          λ              IR                        )                                                                                  ⅆ              log                        ⁢                          xe2x80x83                        ⁢                          I              ⁢                              (                                  λ                  R                                )                                                          ⅆ            t                          ⁢                  (                                                    β                o                            ⁢                              (                                  λ                  IR                                )                                      -                                          β                r                            ⁢                              (                                  λ                  IR                                )                                              )                    -                                                  ⅆ              log                        ⁢                          xe2x80x83                        ⁢                          I              ⁢                              (                                  λ                  IR                                )                                                          ⅆ            t                          ⁢                  (                                                    β                o                            ⁢                              (                                  λ                  R                                )                                      -                                          β                r                            ⁢                              (                                  λ                  R                                )                                              )                    
Note in discrete time                     ⅆ        log            ⁢              xe2x80x83            ⁢              I        ⁢                  (                      λ            ,            t                    )                            ⅆ      t        ≅            log      ⁢              xe2x80x83            ⁢              I        ⁢                  (                      λ            ,                          t              2                                )                      -          log      ⁢              xe2x80x83            ⁢              I        ⁢                  (                      λ            ,                          t              1                                )                    
Using log Axe2x88x92log B=log A/B,                     ⅆ        log            ⁢              xe2x80x83            ⁢              I        ⁡                  (                      λ            ,            t                    )                            ⅆ      t        ≅      log    ⁢          (                        I          ⁢                      (                                          t                2                            ,              λ                        )                                    I          ⁢                      (                                          t                1                            ,              λ                        )                              )      
So, (4) can be rewritten as                                                                         ⅆ                log                            ⁢                              xe2x80x83                            ⁢                              I                ⁡                                  (                                      λ                    R                                    )                                                                    ⅆ              t                                                                          ⅆ                log                            ⁢                              xe2x80x83                            ⁢                              I                ⁡                                  (                                      λ                    IR                                    )                                                                    ⅆ              t                                      ≅                              log            ⁡                          (                                                I                  ⁡                                      (                                                                  t                        1                                            ,                                              λ                        R                                                              )                                                                    I                  ⁡                                      (                                                                  t                        2                                            ,                                              λ                        R                                                              )                                                              )                                            log            ⁡                          (                                                I                  ⁡                                      (                                                                  t                        1                                            ,                                              λ                        IR                                                              )                                                                    I                  ⁡                                      (                                                                  t                        2                                            ,                                              λ                        IR                                                              )                                                              )                                      ≡        R                            (        5        )            
where R represents the xe2x80x9cratio of ratios.xe2x80x9d
Solving (4) for s using (5) gives   s  =                              β          r                ⁢                  (                      λ            R                    )                    -              R        ⁢                  xe2x80x83                ⁢                              β            r                    ⁢                      (                          λ              IR                        )                                              R        ⁢                  (                                                    β                o                            ⁢                              (                                  λ                  IR                                )                                      -                                          β                r                            ⁢                              (                                  λ                  IR                                )                                              )                    -                        β          o                ⁢                  (                      λ            R                    )                    +                        β          r                ⁢                  (                      λ            R                    )                    
From (5) note that R can be calculated using two points (e.g., plethysmograph maximum and minimum), or a family of points. One method using a family of points uses a modified version of (5). Using the relationship                                                         ⅆ              log                        ⁢                          xe2x80x83                        ⁢            I                                ⅆ            t                          =                                            ⅆ              I                        /                          ⅆ              t                                I                                    (        6        )            
now (5) becomes                                                                                                                                     ⅆ                      log                                        ⁢                                          xe2x80x83                                        ⁢                                          I                      ⁡                                              (                                                  λ                          R                                                )                                                                                                  ⅆ                    t                                                                                                              ⅆ                      log                                        ⁢                                          xe2x80x83                                        ⁢                                          I                      ⁡                                              (                                                  λ                          IR                                                )                                                                                                  ⅆ                    t                                                              ≅                              xe2x80x83                            ⁢                                                                                          I                      ⁡                                              (                                                                              t                            2                                                    ,                                                      λ                            R                                                                          )                                                              -                                          I                      ⁡                                              (                                                                              t                            1                                                    ,                                                      λ                            R                                                                          )                                                                                                  I                    ⁡                                          (                                                                        t                          1                                                ,                                                  λ                          R                                                                    )                                                                                                                                  I                      ⁡                                              (                                                                              t                            2                                                    ,                                                      λ                            IR                                                                          )                                                              -                                          I                      ⁡                                              (                                                                              t                            1                                                    ,                                                      λ                            IR                                                                          )                                                                                                  I                    ⁡                                          (                                                                        t                          1                                                ,                                                  λ                          IR                                                                    )                                                                                                                                              =                              xe2x80x83                            ⁢                                                                    [                                                                  I                        ⁡                                                  (                                                                                    t                              2                                                        ,                                                          λ                              R                                                                                )                                                                    -                                              I                        ⁡                                                  (                                                                                    t                              1                                                        ,                                                          λ                              R                                                                                )                                                                                      ]                                    ⁢                                      I                    ⁡                                          (                                                                        t                          1                                                ,                                                  λ                          IR                                                                    )                                                                                                            [                                                                  I                        ⁡                                                  (                                                                                    t                              2                                                        ,                                                          λ                              IR                                                                                )                                                                    -                                              I                        ⁡                                                  (                                                                                    t                              1                                                        ,                                                          λ                              IR                                                                                )                                                                                      ]                                    ⁢                                      I                    ⁡                                          (                                                                        t                          1                                                ,                                                  λ                          R                                                                    )                                                                                                                                              =                              xe2x80x83                            ⁢              R                                                          (        7        )            
Now define
Then describes a cluster of points whose slope of y versus x will give R.
xe2x80x83x(t)=[I(t2,xcexIR)xe2x88x92I(t1,xcexIR)]I(t1,xcexR) y(t)=[I(t2,xcexR)xe2x88x92I(t1,xcexR)]I(t1,xcexIR) y(t)=Rx(t)xe2x80x83xe2x80x83(8)
The optical signal through the tissue can be degraded by both noise and motion artifact. One source of noise is ambient light which reaches the light detector. Another source of noise is electromagnetic coupling from other electronic instruments. Motion of the patient also introduces noise and affects the signal. For example, the contact between the detector and the skin, or the emitter and the skin, can be temporarily disrupted when motion causes either to move away from the skin. In addition, since blood is a fluid, it responds differently than the surrounding tissue to inertial effects, thus resulting in momentary changes in volume at the point to which the oximeter probe is attached.
Motion artifact can degrade a pulse oximetry signal relied upon by a physician, without the physician""s awareness. This is especially true if the monitoring of the patient is remote, the motion is too small to be observed, or the doctor is watching the instrument or other parts of the patient, and not the sensor site.
In one oximeter system described in U.S. Pat. No. 5,025,791, an accelerometer is used to detect motion. When motion is detected, readings influenced by motion are either eliminated or indicated as being corrupted. In a typical oximeter, measurements taken at the peaks and valleys of the blood pulse signal are used to calculate the desired characteristic. Motion can cause a false peak, resulting in a measurement having an inaccurate value and one which is recorded at the wrong time. In U.S. Pat. No. 4,802,486, assigned to Nellcor, the assignee of the present invention, the entire disclosure of which is incorporated herein by reference, an EKG signal is monitored and correlated to the oximeter reading to provide synchronization to limit the effect of noise and motion artifact pulses on the oximeter readings. This reduces the chances of the oximeter locking onto a periodic motion signal. Still other systems, such as the one described in U.S. Pat. No. 5,078,136, assigned to Nellcor, the entire disclosure of which is incorporated herein by reference, use signal processing in an attempt to limit the effect of noise and motion artifact. The ""136 patent, for instance, uses linear interpolation and rate of change techniques to analyze the oximeter signal.
Each of the above-described techniques for compensating for motion artifact has its own limitations and drawbacks. It is therefore desirable that a pulse oximetry system be designed which more effectively and accurately reports blood-oxygen levels during periods of motion.
According to the present invention, a method and apparatus are provided for reducing the effects of motion artifact and noise on a system for measuring physiological parameters, such as, for example, a pulse oximeter. The method and apparatus of the invention take into account the physical limitations on various physiological parameters being monitored when weighting and averaging a series of samples or measurements. Varying weights are assigned different measurements. Optionally, measurements are rejected if unduly corrupt. The averaging period is also adjusted according to the reliability of the measurements. More specifically, a general class of filters is employed in processing the measurements. The filters use mathematical models which describe how the physiological parameters change in time, and how these parameters relate to measurement in a noisy environment; The filters adaptively modify a set of averaging weights and averaging times to optimally estimate the physiological parameters.
In a specific embodiment, the method and apparatus of the present invention are applied to a pulse oximeter which is used to measure the oxygen saturation of hemoglobin in arterial blood. The system takes the natural logarithm of the optical oximetry data and then bandpass filters the data to get absorption-like data. The bandpass filter strongly attenuates data below 0.5 Hz and above 10 Hz in an attempt to remove as much out-of-band noise as possible. This filtered data is then processed through two algorithms: a rate calculator and a saturation calculator.
The system calculates the heart rate of the patient one of three ways using the oximetry data. An adaptive comb filter (ACF) is employed to track the slowly varying heart rate. The tracking of the heart rate by the ACF is quite robust through noisy environments, however, the ACF is not a good heart rate finder. As a result, the system periodically calculates the power spectrum of one of the wavelengths and uses it to find and/or verify the heart rate. In cases of arrhythmia or suddenly changing heart rates, the system employs a pattern matching technique that recognizes sequences of crests and troughs in the data and calculates an average heart rate period over a set number of samples.
The system then employs the calculated heart rate to digitally comb filter the data so that only the energy at integer multiples of the heart rate are allowed through the filter. The comb filter frequency varies as the heart rate varies, attenuating motion energy not at the heart rate or multiples thereof. To remove noise energy at integer multiples of the heart rate, the system adaptively signal averages full cycles of past plethysmographs, i.e., pleths, using a Kalman filter to limit the rate of change in the pleth shape or size.
The system then calculates two saturations, one with the pleth cycle data which has been comb filtered as described above, and one with raw data from the output of the band pass filter. Both saturations are calculated using time based signals and using an adaptive Kalman filter which continuously weights all data according to an estimate of the current noise, and limits the rate of change of saturation to a defined limit (currently 1.3 saturation points per second). Data points that result in a saturation calculation (prior to weighting and averaging) which is obviously not physiologically possible (e.g., negative saturation, or a saturation greater than 100%) are deemed invalid and are not used and are rejected in an xe2x80x9coutliner rejectionxe2x80x9d step in both saturation calculations. The system then arbitrates between the two saturation values based on rules described below to determine the best saturation. For example, the arbitration may be based on such factors as the noise level or the age of the saturation value. The best saturation may also be a weighted average of the different saturation values.
According to a specific embodiment of the invention, a method for reducing noise effects in a system for measuring a physiological parameter is provided. A plurality of measurements is generated corresponding to at least one wavelength of electromagnetic energy transmitted through living tissue. Selected measurements are compared with at least one expected measurement characteristic. A variable weight is assigned to each of the selected measurements based on the comparison, thereby generating a plurality of differently weighted measurements for each wavelength. A first number of weighted measurements is averaged to obtain a filtered measurement, the first number varying according to the manner in which weights are assigned to a plurality of successive weighted measurements. A plurality of filtered measurements are thus generated for each wavelength. The filtered measurements for each wavelength are then combined and calculations resulting therefrom are adaptively filtered using variable weights based on comparing the calculations to an expected calculation. A second number of the weighted calculations are averaged to obtain a filtered calculation, the second number varying according to the manner in which weights are assigned to a plurality of successive weighted calculations.
A further understanding of the nature and advantages of the present invention may be realized by reference to the remaining portions of the specification and the drawings.